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Abstract 

The main purpose of this article is to describe a numerical scheme for 
solving three-dimensional linear Fredholm integral equations of the second 
kind on the cubic domains. The method is based on interpolation by radial 
basis functions (RBFs) based on Gauss-Legendre nodes and weights. Error 
analysis is presented for this method. Finally, several examples are given and 
numerical examples are presented to demonstrate the validity and applica- 
bility of the method. 
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1 Introduction 


Consider the following three-dimensional linear Fredholm integral equation 
of the second kind 


f pad pb 
Heian ’, / / K(e,y, 2%, 8,t)u(r, 9, )drdedt = h(a, y,2),(#,9,2)€D, (1) 
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where h and K are known functions, u(z,y,z) is the unknown function to 
be determined, \ is a constant and D is an cubic domain. 


Integral equations occur in a wide variety of physical applications. They 
are encountered in various fields of science and numerous applications such as: 
elasticity, plasticity, heat and mass transfer, oscillation theory, fluid dynam- 
ics, filtration theory, electrostatics, electrodynamics, game theory, control, 
queuing theory, electrical engineering, economics, medicine, etc. There are 
many different numerical methods for solving integral equations. Computa- 
tional complexity of mathematical operations is the most important obstacle 
for solving integral equations in higher dimensions. The Nystrom method [15] 
and collocation method [5, 14,29] are the most important approaches of the 
numerical solution of these integral equations. 


To avoid the mesh generation, in recent years meshless techniques have 
attracted attention of researchers. In a meshless method, a set of scattered 
nodes are used instead of meshing the domain of the problem [6, 7, 26]. 


Among meshless methods, the radial basis functions (RBFs) method has 
become known as a powerful tool for the scattered data interpolation prob- 
lem. The main advantage of radial basis functions is that they involve a single 
independent variable regardless of the dimension of the problem. One of the 
domain-type meshless methods, the so-called Kansa’s method developed by 
Kansa in 1990 [20,21], is obtained by directly collocating RBFs, particularly 
the multiquadric (MQ), for the numerical approximation of the solution. 
Kansa’s method was recently extended to solve various ordinary and partial 
differential equations including the one-dimensional nonlinear Burgers equa- 
tion [18] with shock wave, shallow water equations for tide and currents sim- 
ulation [17], heat transfer problems [30], and free boundary problems [19,23]. 


Furthermore, the RBFs have been applied on the one-dimensional do- 
mains for solving linear second kind Fredholm and Volterra integral equa- 
tions in [12], linear integro-differential equations in [13], nonlinear Volterra- 
Fredholm-Hammerstein integral equations in [25] and systems of nonlinear 
integral equations in [11]. Also, a numerical solution of two-dimensional Fred- 
holm integral equations of the second kind on the square domains by Gaussian 
radial basis functions without the error analysis is introduced in [1]. 


In this paper, we will use the radial basis functions (RBFs) approxima- 
tion for solving three-dimensional linear Fredholm integral equations of the 
second kind on cubic domains. The remainder of the paper is organized as 
follows: in Section 2, we show that how the radial basis functions are used 
to approximate the solution. In Section 3, we present a numerical method 
for solving the linear Fredholm integral equations of the second kind by the 
RBF approximation. In Section 4, error analysis for the proposed method is 
presented. Numerical examples are given in Section 5. Finally, we conclude 
the article in Section 6. 
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2 An outline of RBFs 


The radial basis function (RBF) method for multivariate approximation is 
one of the most often applied tools in modern approximation theory due to 
spectral accuracy, flexibility with respect to geometry, dimensional indepen- 
dence and ease of implementation especially when the task is to interpolate 
scattered data in multi dimensions. The multiquadric (MQ) method was 
originally introduced by Hardy in 1968 for the interpolation of two dimen- 
sional scattered data to solve a problem from cartography [16]. The problem 
was to construct a continuous function from a set of sparse, scattered mea- 
surements from some source points on a topographic surface, which exactly 
fit to the given data and provides a good approximation of the features of 
the surface such as location of hilltops, saddles, breaks in slope, and drainage 
junctions. In fact the MQ method is a special version of the radial basis func- 
tions method. In 1982 Franke tested a large number of interpolation methods 
for two dimensional scattered data, and found that MQ method was one of 
the most impressive [10]. 


Table 1: Some well-known functions that generate RBFs. 


Name of function Definition 


Multiquadrics (MQ) o(x)=¥ IIxll3 te? 
a 
Inverse multiquadrics (IMQ) o(x) = (y IIx|I3 ar 2) 
1 


Inverse quartics (IQ) o(x) = (Ilxtl2 an eo). 
Gaussian (GA) (x) = exp (—cllxll3) 

Thin plate splines (TPS) (x) =(—1)°" [x|"log [[xIl2 
Conical splines (CS) (x) =[[x||3"* 


Definition 1. [28] A function ¢ : R* —> R is called radial basis provided 
there exists a univariate function ¢ : [0,co) —> R such that 


where r = ||x|| and ||.|| is some norm on R°, usually the Euclidean norm. 


Some well-known RBFs are listed in Table 1. Plots of some RBFs are 
presented in Figure (1). 
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The idea of radial basis function method for interpolation is derived from 
piecewise polynomial interpolation using a function of Euclidean distance and 
defined as follows: 


Figure 1: Plots of some radial basis functions. 


2.1 Radial basis function interpolation 


In the standard RBF interpolation problem, we are given generally scat- 
tered data sites X = {x1,...,xw}C D and associated real function values 
u(x;),72 =1,...,N. Here D is usually some bounded domain in R*. It is our 
goal to find a (continuous) function Y,u : R*° > R that interpolates the 
given data, i.e., such that 


Py,u(x;) =u(x;), t=1,...,N. (2) 


In the RBF literature (see, e.g., [8,28]) one assumes that this interpolant is 
of the form 
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N 
Py u(x) =) Ge(|lx - x5), (3) 


where the basic function y is assumed to be RBF and the coefficients C = 
[c1,...,¢w]" are found by enforcing the interpolation constraints (2). This 
implies that 


AC =U, (4) 


where Ajj = 9(||x; — x;||) and U = [u(xi),...,u(xw)]”. If the ma- 
trix A is such that CTAC is strictly positive for all possible choices of 
X = {x,,...,xw} and C = [cq,...,cn]? € R™ — {0} the solution of the 
interpolation problem is guaranteed. 

The following results establish invertibility of the matrix A for different 
radial basis functions: 


Definition 2. [28] A function ¢ is called completely monotone on (0, co) 
if it satisfies ¢ € C™(0,0o), and 


(-1)'g (r) > 0, 


for all 1 € No and all r > 0. The function ¢ is called completely monotone 
on [0, 00) if it is in addition in C[0, co). 


Theorem 1. /27] If é(r) = p(vr) is completely monotone but not constant 
on [0,00), then for any set of N distinct points {x1,...,xnw}, the Nx N 
matriz A with entries ~(||x; — x,;||) 1s positive definite (and therefore non- 
singular). 


Theorem 2. /24] Let d(r) = y(V/r) € C°[0,00), d(r) > 0 for r > 0, 
a (r) completely monotone but not constant on (0,00), then for any set of N 
distinct points {x1,...,xn}, the N x N matric A with entries ~(||xi — x;||) 
is positive definite. 


It is easy to prove that the radial basis functions IMQ, IQ and GA satisfy 
the sufficient conditions of Theorem 1, whereas the MQ and linear RBFs 
satisfy the sufficient conditions of Theorem 2, and hence for these types of 
RBFs the system (4) is uniquely solvable for any set of distinct data points. 
Although the matrix A is non-singular in the above cases, usually it is very 
ill-conditioned. Therefore, a small perturbation in initial data may produce 
large amount of perturbation in the solution. 

As given in Table 1, the types of RBF are mainly divided into two cat- 
egories, infinitely smooth and piecewise smooth RBFs [2,22]. The infinitely 
smooth RBF's contain a free parameter c, called the shape parameter, which 
affects both the accuracy of a solution and the conditioning of the collocation 
matrix. In Figure 2, a data set is interpolated with the Gaussian function, 
with different shape parameters. A smaller value of c causes the function to 
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become flatter, while increasing c leads to a more peaks RBF. The optimal 
value of the shape parameter that can produce relatively accurate results is to 
be found numerically. But the optimal choice of the shape parameters is an 
open problem which is still under intensive investigation. Several proposals 
for the choice of an adequate shape parameter can be found in the papers of 
Hardy [16], Franke [10] and Fasshauer [9]. All these proposals are somehow 
related with the number of points in the grid and the distance between those 
points. 


Definition 3. Hardy’s shape parameter 
La 
c=0.815d where d= W De di, 
i=1 
where d; is the distance from i*” center to the nearest neighbor and N is the 
total number of centers. 


Definition 4. Franke’s shape parameter 


1.25D 
c= ——, 
VN 


where D is the diameter of smallest circle encompassing all the center loca- 
tions and N is the total number of centers. 


Definition 5. Fasshauer’s shape parameter 
where N is the total number of centers. 


3 Solution of linear integral equations 


In this paper, we solve the three-dimensional linear Fredholm integral equa- 
tion given in the form 


f pad pb 
uen2) rf ff Kenenssthulrs,thdrdsdt = h(e,12),(0.942) € Dy () 
D4 ‘¢. a 


where h and K are known functions, u(z,y,z) is the unknown function to 
be determined, is a constant and D is an cubic domain. 
Consider equation (5) with the following assumptions: 


(i) h £0 is continuous in C(D), 
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Figure 2: The Gaussian function, with different shape parameters, c. 


(ii) K is continuous in C(D x D), 


(iii) kernel K(x, y, z,r,s,t) is real, continuous and bounded in the domain 


D, i.e. 
fo pd pb 
L= _ sup [ f [k@uens.0)drasae < 00. 
(x,y,z)ED e Jc Ja 


Theorem 3. Existence and uniqueness of solution to equation (5) follow by 
assumptions (1) — (itt) and the condition 


|A|Z <1. 


Proof. It can be proved using Banach’s fixed point theorem in a similar 
method as done in [4], Chapter 5 (for one-dimensional linear Fredholm inte- 
gral equations). 


3.1 The proposed method 


To apply the method, we need a RBF ¢ and N nodal scattered points to 
initiate the RBF method. These nodes can be selected arbitrary on the whole 
of the domain D, such as X = {(21, 41, 21),-.-, (en, yn, 2n)}. Therefore, to 
solve equation (5), we estimate the unknown function u(z, y, z) by the RBF 
interpolation method as 
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N 
AT 3 
ula Ye ~ Soe KPK( T,Y,% =C U(x, y, 2), (2, y,z) EDCcR ’ (6) 
k=1 


where 
CT = [G1,---,En], 
Uw" (2, YU; Zz) = [i (x, Y; Z), = ., On (2,Y, z)], 
ox(a,y,2) =9(Ve@ UK)? + (y— ye)? +(z 2), k=1,...,N. 


We replace the expansion (6) with u(x, y,z) and install the collocation 
points (a;, yi, 2:),2 = 1,2,...,N in equation (5). Thus we obtain 


Ppa pe 
oT LW (cities) = rf ‘ / Keiser 8s) (r 8: Odrdsat =h(xi,yi, zi). (7) 
ee iC a 


The integrals in (7) must usually be evaluated numerically. we convert the 
intervals [a, b], [c,d] and [e, f] to the interval [—1, 1] by using a simple linear 
transformations of the form 


b-a b+a b-a 
SEs ae 5 = 9(f) => dr = ——df, 
d—c d+e d—c 

$= > n+ —— = h(n) = ds= 5 
f-e fte f-e 
t= Ber 5 m(r) dt 5 dr, 


and so equation (7) takes the following form: 


on fy Ci wf is as K (wi, vir 21, 9(8), b(n), m(7))W (a(6),Ma)som(r) ganar} 


= A(x, Yas 2)s 
(8) 


where 


_ A(b—a)(d—o)(f —e) 
= 7 


Using an my -point Gauss - Legendre quadrature formula with the points 
Tp, 8q, tx in the interval [—-1,1] and weights w,, wg, wx for numerical integra- 
tion in equation (8), we can approximate the integral 


I [ K(a2,y, 2, 9(€), h(n), m(7)) U(g(€), h(n), m(7) )dEdndr, (9) 


with 


Radial basis functions method for solving three-dimensional ... 23 


mn MN MN 


WpWqwek (a, Y,4, gE), hq); m(TK))¥(9(Ep), h(nq); m(Tk)). 


(10) 


Utilizing this numerical integration rule in equation (8), we obtain the 
following linear system of algebraic equations 


Kueeds ais 
oF fw(es, wer 2) HDS DoD wre ttak C2 tes 21 906): Cg) 2m(74)) HOLE): MC) m(re))} 


= h(xi, yi, zi), 
(11) 
where i = 1,2,...,N. This is a linear system of equations that can be solved 
by iterative methods to obtain the unknown vector C7. 


4 Error analysis 


This section includes the error estimate and the rate of convergence of the 
presented method. To understand the numerical behavior of the interpolant 
or approximant it is essential to have bounds on the approximation error 
and on the condition number of the interpolation matrix. These bounds 
are usually expressed employing two different geometric measures. For the 
approximation error, it is crucial to know how well the data sites X fill the 
region . This can be measured by the fill distance 


Psy 5 = ae Pesta I|x = xjll2, 
which gives the radius of the largest data-site free ball in D. The condition 
number, however, will obviously only depend on the data sites X and not on 
the region D. Moreover, if two data sites tend to coalesce then the correspond- 
ing interpolation matrix has two rows which are almost identical. Hence, it 
is reasonable to measure the condition number in terms of the separation 
distance 


= min || | 
— in ||x; — x,]]. 
dx 2 ixj - 


A set X of data sites is said to be quasi-uniform with respect to a constant 
Cy, > 0 if 


qx S hyp < Coudx: 


Definition 6. [28] The definition of the native space is 
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Ng(D) = {! € Io(R*)NC(R*): tee rae} 


vé 
where ¢ is a Fourier transform of ¢. 


Theorem 4. /2/ Let ¢ is positive definite RBF with infinitely smoothness. 
Suppose that D C R® be open and bounded, satisfying an interior cone con- 
dition. Denote the interpolant of a function u € X¢(D) based on this RBF 
and the distinct set X = {x,...,an} by P,u. Then for every 1 € N there 
exist constants ho(l),C; such that 


|lu — Py ulln~(p) < Crh), , llullys(v), (12) 


for all x € D, provided hy, < ho(I). 


Remark 1. As a conclusion from Theorem 4, for Gaussians ¢(x) = 
e(—ellxll?) c > 0, we get for some positive constant | that 


(ee) 
RxD | 


|u — Py ullz-o(p) < e |ull4(D)> (13) 


provided that h, ,, is sufficiently small and u € Ng(D). 
The corresponding result for (inverse) multiquadrics ¢(x) = (||x|| + 
c*)*,c>0,a <0, ora>0anda4¥N, is 


Paes ee 
lu — Pyullz-(p) < elm ) Ilellx4(D)> (14) 
For thin plate splines ¢ (x) =(—1)**"||x||”“log |x|, , & € N, we get 
Iu — Pyull ssc) $ CRE, llulhsgco) (15) 


Let #% be the Urysohn integral operator: 


f pd pb 
(Huylay2)=A i: i: K(a,y, 2,7, 8, t)u(r, s, t)drdsdt, (16) 


we can rewrite the integral equation (5) in operator form as 
u- #Hu=h. (17) 


Define the approximating operator #y,N > 1, on C(D) by 


mn TN TAN 


Anu, Y z) =- oe > > wpwqurk (x, Y,2, g(Ep); h(nq); m(tk))C7 U(g(Ep), h(nq); m(TK)). 
p=l1lq=1k=1 
(18) 


The abstract form of equation (8) is 


Radial basis functions method for solving three-dimensional ... 25 
un — Py, Hun = Ah, 
for N sufficiently large and the above equation can be rewritten as 


(l- P, Huy = Pyh. (19) 


Remark 2. Note that, Y, : C(D) — Vy is the collocation projection 
operator on the collocation points X = {(21,y1, 21),--.,(@n, yn, zn) } C D, 
where the subspace Vy := span{¢1,...¢n} C C(D) has finite dimension. 


Suppose t € Vx, we simply have Y,, u. 


a = 
The abstract form of equation (11) is 


which shows that the scheme is a discrete collocation method [3]. 
Consequently an iterated discrete collocation solution can be obtained. For 
this purpose we set 


tn =h+ 2, An (tn), (21) 


and by applying the operator Y,, on both sides of (21), and using the relation 
(20) we simply have 


P tin = tin. (22) 
Thus we conclude 
(I -— Py, Hn )Un =h. (23) 


Theorem 5. /3] Assume the family {4} of (8) is collectively compact 
and pointwise convergent on C(D). Let {Y,,} be a family of interpolatory 
projection operators on C(D) to C(D), and assume 


Pyu—u as N—-+o, (24) 


for allu€ C(D). Finally, assume the integral equation (17) is uniquely solv- 
able for allh € C(D) and ux be a unique solution of this equation. Then for 
all sufficiently large N, say N > M, (I— #y P,,)~! exists and is uniformly 
bounded. Also, for the solution tn 

lin — Ulla (Dy SI - An Py) MU — Av Pytielln~(p)- (25) 


Remark 3. From (24) and the principle of uniform boundedness for the 
radial basis functions, 


cy = sup | Py || < 00. (26) 
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Similarly, from the pointwise convergence of { #y}, 


Ck = sup ||. Z|] < 00. (27) 


Theorem 6. Having in mind the assumptions of (5). Suppose that us € 
Ng(D) is the unique exact solution of equation (5) and the proposed method 
has been installed on the quasi-uniform set X = {(xi, yi, 2) }§_,. Then there 
exists M > 0 such that for every N > M, the method has a unique solution 
UN 


|| — us|[n~e(p) +0 as N —+ oo. (28) 
In addition, for the iterative solution tin for equation (23) we have 
un — Uelln~(D) < er{|| Hux — Hn usl|b-(p) + ex(1 + Cp) Crh", || Uelle.,(D)}> 


provided that u, € X¢(D), and for the discrete collocation solution tin of 
equation (20) we have 


Gn — uUxllnc0(D) S cper|| 4% ux — Hn ux||peocp) + 1 + ep) + eperce Crh. [lulls (D)> 


where cy < co is a bound for (I -— #y Py) 71. 


Proof. From Theorem 5, the iterated method has a solution ty and 
ln — wxllze(D) < IE -— Hn Py) "I Xue — HN Py Uxll 20d) 
<c7||Hus. — #n Px Ux|| L200 (dD) 
< ¢7 {|| ux — Hy ux||re0(p) + || #n (ux — Py Ux)||b20 (pd) } 
<e7{|| Hux — Hn ux|| boca) + ckl| (ux — Ay Ux) ll be (dy } 
< cr{|| Hux — An ux||p-0(p) + ck (lux — G«|| e0(D) + || Py Ue — Py G«|| E00 (D))} 


<er{|| Hux — Hn us||pocp) + ce(1 + ep)||(ux — @)|| 200 (Dy } 


<er{|| us — #nus||pocp) + ee(1 ep)Crh\, 5 Ilus lle, (D)} 


The last inequality is implied by (15). Moreover, let in = pntin, and 
consider the decomposition 


Ux — in = Us — Pytin = (Ux — Ayus) + Ay (ux — tin), (29) 


which yields 
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I[fiv — Us I Z2°(D) S ||Ue — Pry Ul] r2¢(D) + plan — Us|] 2-0 (D) 
<(1+ Cp) Crh’, [ele (D) + p(cr{|| Hux — An us||b~(d) 
+cn(1+ Cp) CIR’, |ux|4(D) }) 
< cpez|| Hux — Hy ux||_ t(D) 
So See 


+ (1+ ep)(1 + cpercr) Crh), , [lus lle4(D) 
_-_-_-.uaaau“rew——— 


N—+o0 hy p— 0 


Finally, we obtain 


|| — us|| t(D) —> 0 as N—>+o. 


Corollary 1. Theorem 6 shows that, both the quadrature and the RBF 
approximation error bounds affect the final estimation. If for a sufficiently 
smooth kernel K(ax,t,y,s) a high order quadrature is employed then the total 
error is dominated by the error of the RBF approximation. 


5 Numerical examples 


In this section, we present some numerical examples where D is a bounded 
domain in R*. In addition to the Hardy, Franke and Fasshauer shape pa- 
rameters, the minimum error obtained by trial and error is presented. We 
have used the ten-point Gauss-Legendre quadrature rule for numerical inte- 
gration. All of the computations have been done using the Maple 14 with 
just 80 digits precision. we calculate the RMS error in 2197 points that are 
distributed uniformly in the computational domain. The RMS error of the 
numerical result is described using 


1 2197 
= =e ss . eee : Fs =) 2 
RMS = 2197 2 | Ui Deets 2s) Ol Diba, 26) | 


where u(x, y, z) is the exact solution, &(2, y, z) is the approximate solution. 


Example 1. Consider the three-dimensional linear Fredholm integral equa- 
tion 


rep pl ayz 
u(x, Yy,z)— = 4 — u(r, s,t)drdsdt = f(x,y, z), 
@ua—3f ff pe curso flay) 


where 
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LYZ 
2+ 27 + 2y +22 


f(x,y, 2) = sin? (x) sin?(y) sin?(z) + (cos(1) sin(1) 1)° 


The exact solution for this equation is u(2, y, 2) = sin?(a) sin?(y) sin?(z). 


Chebyshev 


Figure 3: Node distribution with 216 nodes for Example 1. 


The results obtained by the Chebyshev distribution of points for different 
numbers of N in terms of RMS is given in Tables 2, 3 and 4. The distribution 
of nodes are depicted for N = 216 in Figure 3. As we expected, from Theorem 
6, the results converge to the exact values along with the increase of the nodes. 
In computations, for GA RBFs, to obtain the better result, we can use the 
small(big) parameter c but the condition number of the final system is grown 
fast instead. As can be seen, the convergence rate of the method is arbitrarily 
high when GAs and IQs are used and is approximately Oe _») when TPSs 
and CSs are used. Therefore the numerical results confirm the theoretical 
error estimates. 

Example 2. Consider the three-dimensional linear Fredholm integral equa- 
tion 


2 pl pl 
u(x, y, Zz) — aL | [) siney2)ulr, s,t)ardsdt = f(z,y, 2), 


where 


f(a,y,2) = 22? — = sin(2y2), 


with exact solution u(x, y, 7) = x?y?z?. The results obtained by the Regular 
distribution of points for different numbers of N in terms of RMS is given 
in Tables 5, 6 and 7. The distribution of nodes are depicted for N = 294 
in Figure 4. It should be noted that for a fixed sufficiently large my, by 
increasing N, the error of the method is of O(n.) (i.e., arbitrarily high 


for the use of MQs and IMQs and approximately OnE 5) for the use of 
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Table 2: Results provided by Chebyshev distribution for Example 1. 
N Bigs ds shape Cc GA 
parameter 
Hardy 0.70 786.5610" 
27 0.7339 Franke 0.72 8.13 K 1057 
Fasshauer 0.38 2.52 x 1074 
a* 0.38 2.52 x 1074 
Hardy 0.44 7.88 x 1074 
64 0.6334 Franke 0.50 1.59 x 1073 
Fasshauer 0.25 4.83 x 10-3 
a* 0.43 Ta Os 
Hardy 0.30 2.11 x 1074 
125 0.4991 Franke 0.37 2.51 x 1074 
Fasshauer 0.18 7.50 x 107° 
a* 0.20 ADA sc 10-° 
Hardy 0.22 5.13 x 107° 
216 0.4189 Franke 0.28 2.02 x 1074 
Fasshauer 0.14 2.35 x 1074 
a* 0.20 2.20 x 1076 
a* indicates the shape parameter which minimum error is observed. 
Table 3: Results provided by Chebyshev distribution for Example 1. 
WO hee Cs TPS 
k=1 k=2 k=1 k=2 
27 0.7339 9.14 x 107? 3.92107? 1.56107! 2.21 x 107! 
64 0.6334 3.36 x 107? 1.71 x 107? 5.88x 107? 4.82 x 107? 
125 0.4991 2.29 x 107? 5.76 x 1073 1.17x 107? 4.54 x 1073 
216 0.4189 1.71 x 107? 3.34 x 107% 6.05 107% 2.64 x 107-3 
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Table 4: Results provided by Chebyshev distribution for Example 1. 
N Baad shape c IQ 
parameter 
Hardy 0.70 2.24 x 107? 
27 0.7339 Franke 0.72 2.29 x 107 
Fasshauer 0.38 3.52 x 1072 
a* 1.25 1.57 x 1073 
Hardy 0.44 7.44 x 1073 
64 0.6334 Franke 0.50 6.27 x 1073 
Fasshauer 0.25 9.30 x 1073 
a* 0.93 8.91 x 10-4 
Hardy 0.30 6.60 x 1073 
125 0.4991 Franke 0.37 4.82 x 1073 
Fasshauer 0.18 8.22 x 1073 
a* 1.61 1.43 x 1074 
Hardy 0.22 4.34 x 1073 
216 0.4189 Franke 0.28 2.35 x 1073 
Fasshauer 0.14 6.23 x 1073 
a* 0.92 7.84 x 107° 


a* indicates the shape parameter which minimum error is observed. 


TPSs and CSs ) because the RBF interpolation error overcomes the error of 
integration method and so increasing N has no significant effect on the error. 
Based on the obtained fill distance, the results confirm Theorem 6. 


Example 3. Consider the three-dimensional linear Fredholm integral equa- 
tion 


tpt ph 
u(ays2)— 5 f | i xyzu(r, s,t)drdsdt = f(x,y, 2), 


where 


1 soo 
f(z, y, 2) = xyz exp( a y” z*) 16 (1 =) TYZ, 
with exact solution u(z, y, z) = xyz exp(—a? — y? — 27). The results obtained 
by the Halton distribution of points for different numbers of N in terms of 
RMS is given in Tables 8, 9 and 10. A set of 216 Halton points in the unit 
cubic in R® has been shown in Figure 5. For sufficiently large my, increasing 
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Figure 4: Node distribution with 294 nodes for Example 2. 
Table 5: Results provided by Regular distribution for Example 2. 
N bh iy shape Cc MQ 
parameter 
Hardy 0.41 1.01 x 10-? 
48 0.5246 Franke 0.54 9.82 x 1073 
Fasshauer 0.29 1.12 x 107? 
ax 1.82 7.21 x 1073 
Hardy 0.27 5.16 x 1073 
100 0.3842 Franke 0.37 4.81 x 107% 
Fasshauer 0.20 5.26 x 1073 
a* 4.23 9.39 x 1074 
Hardy 0.20 3.26 x 107% 
180 0.3039 Franke 0.28 3.09 x 1073 
Fasshauer 0.15 3.30% 10-3 
a* 5.04 2.20 x 1074 
Hardy 0.16 2.16 x 1073 
294 0.2517 Franke 0.22 2.04 x 1073 
Fasshauer 0.12 2.19 x 1073 
a* 5.52 3.88 x 107° 


a* indicates the shape parameter which minimum error is observed. 


the number of integration nodes my has no significant effect on the error and 
the proposed method will be of O(h). .), by increasing N. The numerical 
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Table 6: Results provided by Regular distribution for Example 2. 
N 3 shape Cc IMQ 
parameter 
Hardy 0.41 1.35 x 107? 
48 0.5246 Franke 0.54 1.39 x 107? 
Fasshauer 0.29 1.20 x 107? 
a* 2.79 4.72 x 1073 
Hardy 0.27 7.23 x 1073 
100 0.3842 Franke 0.37 7.73 x 1073 
Fasshauer 0.20 6.32 x 1073 
a* 4.71 6.68 x 1074 
Hardy 0.20 4.51 x 1073 
180 0.3039 Franke 0.28 5.08 x 1073 
Fasshauer 0.15 4.11 x 1073 
a* 8.40 8.02 x 10~° 
Hardy 0.16 3.07 x 1073 
294 0.2517 Franke 0.22 3.51 x 1073 
Fasshauer 0.12 3.13 x 1073 
a* 5.52 3.88 x 107° 


a* indicates the shape parameter which minimum error is observed. 


Table 7: Results provided by Regular distribution for Example 2. 


‘ne = Cs TPS 
k=1 k=2 k=1 k=2 
48 0.5246 9.92 x 10-2 8.70 x 10-2 1.05x 107! 1.00 x 1074 
100 0.3842 4.87 x 1072 3.76 x 10-2 4.63x 1072 3.15 x 107? 
180 0.3039 3.18 x 1072 2.06 x 10-2 2.71 x 10-2 1.53 x 107? 
294 0.2517 92.22 x 1072 1.16 x 10-2 1.69x 10-2 7.76 x 1073 
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results confirm the theoretical error estimates. In addition, the numerical 
results shows that the accuracy of RBFs with shape parameter is better than 
RBFs without shape parameter. 


Figure 5: Node distribution with 216 nodes for Example 3. 


6 Conclusion 


In this paper, a collocation method based on RBFs for numerical solution of 
three-dimensional linear Fredholm integral equations of the second kind on 
the cubic domains is presented. The proposed method is a meshless method, 
which requires no domain elements for the interpolation or approximation. 
The Gauss-Legendre quadrature formula is employed for numerical integra- 
tion. Error analysis was provided for suffficiently smooth kernel and source 
functions. The method is very convenient for solving higher dimensional in- 
tegral equations because the RBF is defined as the function of distance. The 
proposed method can be easily expanded into non-cube domains. It means 
the process of solving is no more complicated in spite of increasing the dimen- 
sion of problem. This is significant advantage of the method in comparision 
with other strategy for solving Integral Equations in three dimension. 
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Table 8: Results provided by Halton distribution for Example 3. 


N Big ds shape Cc GA 
parameter 
Hardy 0.22 4.16 x 1073 
27 0.4912 Franke 0.29 3.88 x 1073 
Fasshauer 0.38 3.56 x 1073 
a* 1.57 1.81 x 1073 
Hardy 0.14 9.10 x 10-4 
64 0.3456 Franke 0.20 7.84 x 1074 
Fasshauer 0.25 6.89 x 1074 
a* 0.88 8.57 x 107° 
Hardy 0.12 1.76 x 107-4 
125 0.2820 Franke 0.17 1.51 x 1074 
Fasshauer 0.18 1.46 x 10-4 
a* 0.84 VAD «107° 
Hardy 0.09 2.14 x 1075 
216 0.2780 Franke 0.13 1.70 x 1075 
Fasshauer 0.14 1.61 x 10-5 
a* 0.90 9.10 x 1078 


a* indicates the shape parameter which minimum error is observed. 


Table 9: Results provided by Halton distribution for Example 3. 


‘ne = Cs TPS 
a k=1 k=2 k=1 k=2 
27 (0.4912 5.34 x 1072 1.09 x 10-2. 2.12x 10-2 4.07 x 107-2 
64 0.3456 3.61 x 1073 3.22x 10-3 147x107? 5.90 x 10-3 
125 0.2820 2.28 x 1073 1.60 x 10-3 9.151073 3.43 x 1073 
216 0.2780 1.74 x 1073 4.88x10-4 2.041073 1.35 x 1073 
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Table 10: Results provided by Halton distribution for Example 3. 
N eee shape c MQ 
parameter 
Hardy 0.22 3.91 x 1073 
27 0.4912 Franke 0.29 3.77 x 1073 
Fasshauer 0.38 3.65 x 1073 
a* 1.60 2.95 x 1073 
Hardy 0.14 2.38103 
64 0.3456 Franke 0.20 2.14 x 1073 
Fasshauer 0.25 1.99 x 1073 
a* 1.80 2.06 x 10-4 
Hardy 0.12 164-x.10-3 
125 0.2820 Franke 0.17 1.46 x 1073 
Fasshauer 0.18 1.43 x 1073 
a* 2.60 1.74 x 107° 
Hardy 0.09 1.20 x 1073 
216 0.2780 Franke 0.13 1.02 x 1073 
Fasshauer 0.14 9.73 x 1074 
a* 2.10 1.46 x 1076 


a* indicates the shape parameter which minimum error is observed. 
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